% %%% SC 07 %%%%


\section {Octree Meshing}
\label{sec:oda}

By octree meshing we refer to the construction of a data structure on top of the linear octree that allows FEM type calculations. In this section, we describe
how we construct the support data structures in order to perform the matrix-vector multiplications  (\texttt{MatVec}) efficiently. The data structure is designed
to be cache efficient by using a Morton ordering based element traversal, and by reducing the memory footprint using compressed representations for both the
octree and the element-node connectivity tables. The algorithm for generating the mesh given a distributed, sorted, complete, balanced linear octree is outlined
in Algorithm \ref{alg:meshing}.

\begin{table} 
\centering
\rule{\textwidth}{0.01mm}
\begin{algorithm}{ \textsc{Octree Meshing And Compression}}\\ 
\rule{\textwidth}{0.01mm}
\flushleft
{\bf Input:~} A distributed sorted complete balanced linear octree, $L$\\ 
{\bf Output:} Compressed Octree Mesh and Compressed Octree.\\
~\\
1. Embed $L$ into a larger octree, $O$, and add boundary octants. \\
2. Identify `hanging' nodes.\\
3. Exchange `Ghost' octants.\\
4. Build lookup tables for first layer of octants. (Section \ref{sec:exhaustive}) \\
5. Perform 4-way searches for remaining octants.  \\
~~~~(Section \ref{sec:negative})\\
6. Store the levels of the octants and discard the anchors.\\
7. Compress the mesh (Section \ref{sec:compress}).
\label{alg:meshing}
\end{algorithm}
\rule{\textwidth}{0.01mm}
\end{table}

\subsection{Computing the Element to Node Mapping}
\label{sec:buildLut}
The nodes (vertices) of a given element are numbered
according to the Morton ordering. An example is shown in Figure
\ref{fig:nodes}. The 0-node of an element, the node with the
smallest index in the Morton ordering, is also referred to as the {\em
``anchor''} of the element. An octant's configuration with respect to
its parent is specified by specifying the node that it shares with its
parent. Therefore, a 3-child octant is the child that shares its third node
with its parent. Nodes that exist at the center of a face of another
octant are called face-hanging nodes. Nodes that are located at the center of 
an edge of another octant are called edge-hanging nodes.

Since all nodes, except for boundary nodes, can be uniquely associated with an element (the element with its anchor at the same coordinate as the node) we use an interleaved representation where a common index is used for both the elements and the nodes. Because of this mapping, the input balanced octree does not have any elements corresponding to the positive boundary nodes. To account for this we embed the input octree in a larger octree with maximum depth $D_{max}+1$, where $D_{max}$ is the maximum depth of the input octree. All elements on the positive boundaries in the input octree add a layer of octants, with a single linear pass ($\mathcal{O}( n / p )$), and a parallel sort ($\mathcal{O}( n / p\log n / p )$)).
The second step in the the computation of the element-to-node mapping
is the identification of hanging nodes.  All nodes which are not
hanging are flagged as being {\em nodes}. 
Octants that are the $0$ or $7$ children of their parent
($a_0, a_7$) can never be hanging; by default we mark them as {\em
nodes} (see Figure \ref{fig:nodes}). Octants that are $3,5,6$ children of their parent ($a_3, a_5,
a_6$) can only be {\em face hanging}, and their status is determined
by a single negative search.\footnote{By {\em ``negative''} searches
we refer to searches in the $-x$ or $-y$ or $-z$ directions. We use
{\em ``positive searches''} to refer to searches along the positive
directions.} The remaining octants ($1,2,4$ children) are {\em edge
hanging} and identifying their status requires three searches.

After identifying hanging nodes, we repartition
 the octree using the algorithm described in Section \ref{sec:blkPart} and all octants
  touching the inter-processor boundaries are communicated to the neighbouring processors. On the processors
   that receive such octants, these are the ghost elements, and their corresponding nodes are called ghost nodes. In our implementation
   of the \texttt{MatVec}, we do not loop over ghost elements recieved from a processor with greater rank and we do not write to ghost values. This framework gives rise to a subtle special case for \textit{singular blocks}. A singular block is a block (output of the partition algorithm), which is also a leaf in the underlying fine octree (input to the partition algorithm). If the singular block's anchor is hanging, it might point to a ghost node and if so this ghost node will be owned by a processor with lower rank. This ghost node will be the anchor of the singular block's parent. We tackle this case while partitioning by ensuring that any singular block with a hanging anchor is sent to the processor to which the first child of the singular block's parent is sent. We also send any octant that lies between (in the Morton ordering) the singular block and its parent to the same processor in order to ensure that the relative ordering of the octants is preserved.
   
After exchanging ghosts, we perform independent sequential searches on each processor to build the element-to-node mappings. We present two methods for the same, an exhaustive approach (Section \ref{sec:exhaustive}) that searches for all the 8 nodes for each element, and a more efficient approach that utilizes the mapping of its negative face neighbours (Section \ref{sec:negative}) to construct its own mapping.
   
\subsubsection{Exhaustive searches to compute mapping}
\label{sec:exhaustive}

The simplest approach to compute the element-to-node mapping would be to
search for the nodes explicitly using a parallel search algorithm, followed by 
the computation of global-to-local mappings that are necessary to manage the distributed data. 
However, this would incur expensive communication and synchronization costs.
To reduce these costs, we chose to use {\em a priori} communication of ``ghost'' octants\footnote{The use of blocks makes it easy to identify ``ghost'' octants.} 
followed by independent local searches on each processor that require no communication. A few 
subtle cases that can not be identified easily during a priori communication are identified
 during the local searches and corrected for later. 

For any element, all its nodes are in the positive direction (except the anchor).
The exhaustive search strategy is as follows: Generate search keys at the location of
 the eight nodes of the element at the maximum depth and search for them in the linear
  octree. Since the linear octree is sorted, the search can be performed in $\mathcal{O}(\log n)$. If the search result is a {\em node}, then the lookup table is updated. If we discover a hanging node instead, then a secondary search is performed to recover the correct (non-hanging) node index. As can be seen from Figure \ref{fig:nodes}, hanging nodes are always mapped to the corresponding nodes of their parent\footnote{The 2:1 balance constraint ensures that the nodes of the parent can never be hanging.}. Unfortunately, secondary searches can not be avoided despite the identification of hanging nodes prior to searching. This is because only the element whose anchor is hanging knows this information and the other elements that share this vertex must first search and find this element in order to learn this information.

\begin{figure*}
  \begin{center} \subfloat[] { \label{fig:nodes}
  \includegraphics{figures/hanging} } \hspace{0.1in} \subfloat[] {
  \label{fig:ghosts} \includegraphics{figures/ghosts} } 
  \caption{{\texttt (a)} Illustration of nodal connectivities required to perform conforming
  finite element calculations using a single tree traversal. Every octant has at least 2 non-hanging nodes, one of which is
  shared with the parent and the other is shared amongst all the
  siblings. The octant shown in blue ($a$) is a child 0, since it shares its zero node
  ($a_0$) with its parent. It shares node $a_7$ with its siblings. All other nodes, if hanging, 
  point to the corresponding node of the parent octant instead. Nodes, $a_3, a_5, a_6$ are face
  hanging nodes and point to $p_3, p_5, p_6$, respectively. Similarly
  $a_1, a_2, a_4$ are edge hanging nodes and point to $p_1, p_2,
  p_4$. {\texttt (b)} The figure explains the special case that occurs
  during exhaustive searches of ghost elements. Element anchored at $a$,
  when searching for node $b$, will not find any node. Instead, one of
  the hanging siblings of $b$, ($c$, $d$) which are hanging will be
  pointed to. Since hanging nodes do not carry any nodal information,
  the nodal information for $b$ will be replicated to all its hanging
  siblings during update for the ghosts.  } 
  \end {center}
\end{figure*}

Using exhaustive searches we can get the mapping for most elements, but
certain special cases arise for ghost elements. This case is illustrated in Figure \ref{fig:ghosts}, where
the ghost elements are drawn in red and the local elements are drawn
in blue. Consider searching for the $+z$ neighbor of element
$a$. Since $a$ is a ghost, and we only communicate a single layer of
ghost octants across processor boundaries, the node $b$ will not be
found. In such cases, we set the mapping to point to one of the
hanging siblings of the missing node ($c$ or $d$ in this case). The
most likely condition under which $b$ in not found is when the $+z$
neighbor(s) is smaller. In this case, we know that at least one of
the siblings will be hanging. Although the lookup table for this
element is incorrect, we make use of the fact that the lookup table
points to a non-existent node, and the owner of the true node simply
copies its data value to the hanging node locations, thereby ensuring
that the correct value is read. This case can only happen for ghost
elements and need not be done for local elements.  In addition, we
occasionally observe cases where neither the searched node nor any of
 its siblings is found.  Such cases are flagged and at the end of the lookup
  table construction, a parallel search is done to obtain the missing nodes
   directly from the processors that own them.

\subsubsection{Four-way searches to compute mapping}
\label{sec:negative}
The exhaustive search explicitly searches for all nodes and in many
cases is the only way to get the correct element to node
mapping. However, it requires a minimum of 7 and a maximum of 13 searches per element. All searches are performed on a sorted list, and can
be done in $\mathcal{O}(\log n)$. In order to reduce the constants
associated with the exhaustive search, we use the exhaustive search only for the
 first layer of octants (octants that do not have neighbours in the negative $x,y$ and $z$ directions). 
For all other octants, the lookup 
information can be copied from the elements in the negative directions. 
Each element in the negative $x,y$ and $z$ directions that shares a face 
with the current element, also shares 4 nodes. Therefore, by performing negative
searches along these directions, we can obtain the lookup information
for 7 out of the 8 nodes of an element. Only the last node, labeled
$a_7$ in Figure \ref{fig:nodes}, cannot be obtained using a negative
search and a positive search is required.

\begin{figure*}
  \begin{center}
	\subfloat[] {
		\label{fig:newSearch}
		\includegraphics{figures/negSearch} }
	\hspace{0.1in}
	\subfloat[] {
		\label{fig:newSwap}
		\includegraphics{figures/negSwap} }		
	\caption{Computing element to node mapping using negative searches. 
	{\texttt (a)} If the found octant ($a$) is not a sibling of the current octant ($b$),
	then the element to node mapping can be copied via the mapping $b_0\leftarrow a_2$, $b_1\leftarrow a_3$,
	$b_4\leftarrow a_6$, and $b_5\leftarrow a_7$. {\texttt (b)} In case the
	found octant ($a$) is a sibling of the current octant ($b$), then the
	mapping depends on whether or not the node in question is
	hanging. If the node is not hanging, then the same mapping as used
	in {\texttt (a)} can be applied. If the node is hanging, then the
	corresponding node indices for the found element are directly
	copied. For the case shown, $(b_0,b_2,b_4,b_6)\leftarrow (a_0, a_2,
	a_4, a_7) = (p_0, p_2, p_4, a_7)$. }  
	\end {center}
\end{figure*}

In order to get the mapping information using negative searches, we
perform the search in the negative direction and check if the current element
is a sibling of the element obtained via the negative search. If the element found by the search
 is not a sibling of the current element, then the lookup information can be copied via a
mapping. For the example shown in Figure \ref{fig:newSearch}, given the
element $b$ and searching in the $-y$ direction, we find $a$, then the node
mapping is $(b_0,b_1,b_4,b_5)$ = $(a_2,a_3,a_6,a_7)$. Corresponding
mappings are $(b_0,b_2,b_4,b_6)$ = $(a_1,a_3,a_5,a_7)$, and
$(b_0,b_1,b_2,b_3)$ = $(a_4,a_5,a_6,a_7)$, for negative searches along
the $x$ and $z$ axes, respectively. Unfortunately, the mapping is a bit more complex 
if the negative search returns a sibling of the current element. If the node in
question is not hanging, then we can copy its value according to the
above mentioned mapping. However, if the node in question is hanging,
then instead of the mapping, the corresponding node indices from
element $a$ are copied. This case is explained in Figure
\ref{fig:newSwap}, where we observe that if node $a_1,b_0$ is hanging,
we need to use $b_0=a_0$ and use $b_0=a_1$ if it is a node.

%Since we have global information about whether a node is hanging or not, this
%swap can be performed easily.

\subsection {Mesh Compression}
\label{sec:compress}
One of the major problems with unstructured meshes is the storage overhead. In the case of the octree, this
 amounts to having to store both the octree and the lookup table. In order to
reduce the storage costs associated with the octree, we compress both
the octree and the lookup table. The sorted, unique, linear
octree can be easily compressed by retaining only the offset of the first
element and the level of subsequent octants. Storing the offset for
each octant requires a storage of three integers (12 {\texttt bytes}) and
a {\texttt byte} for storing the level. Storing only the level represents
a 12{\texttt x} compression as opposed to storing the offset for every
octant.

It is much harder to compress the element to node mapping, which requires eight integers for each element. 
In order to devise the best compression scheme, we first estimate
the distribution of the node indices. The following lemma helps us
analyze the distribution of the indices of the nodes of a given
element.

\begin{lemma}
  The Morton ids of the nodes of a given element are greater than or equal to
  the Morton id of the element.
\end{lemma}
\begin{proof}
  Let the anchor of the given element be $( x, y, z )$ and let its size be
  $h$. In that case the anchors of the 8 nodes of the element are given by $(
  x, y, z ), ( x + h, y, z ), ( x, y + h, z ), ( x + h, y + h, z ) \cdots$ .
  By the definition of the Morton ordering all of these except $( x, y, z )$ are
  greater than the Morton id of the element. The node at $( x, y, z )$ is
  equal to the Morton id of the element.
\end{proof}

\begin{corollary}
 Given a sorted list of Morton ids corresponding to the combined list
 of elements and nodes of a balanced linear octree, the indices of the
 8 nodes of a given element in this list are strictly greater than the
 index of the element. Moreover, if the nodes are listed in the Morton
 order, the list of indices is monotonically increasing. If we
 store offsets in the sorted list, then these offsets are strictly
 positive.
\end{corollary}

%\begin{corollary}
%  If the indices are stored as offsets from the index of the given element,
%  then the offsets are not sensitive to the distribution of octants whose
%  Morton id is less than that of the given element.
%\end{corollary}

Based on these observations we can estimate the expected range of offsets. 
Let us assume a certain balanced
octree, $O$, with $n$ octants (elements and hanging-nodes) and with 
maximum possible depth $D_{\max}$. Consider an element in the
octree, $o_i$, whose index is $i$, $0 \leqslant i < n$. The offset of
the anchor of this element is either $i$ (if the anchor is not hanging)
or $n_0<i$.  The indices for the remaining 7 nodes do not depend on 
octants with index less than $i$. In addition
since the indices of the 7 nodes are monotonically increasing, we can
store offsets between two consecutive nodes. That is, if the indices
of the 8 nodes of an element, $o_i$, are $( n_0, n_1, n_2, n_3, n_4,
n_5, n_6, n_7 )$, we only need to store $( n_0-i, n_1 - n_0, n_2 -
n_1, n_3 - n_2, n_4 - n_3, n_5 - n_4, n_6 - n_5, n_7 - n_6 )$. 
To efficiently store these offsets, we need to estimate how
large these offsets can be. We start with a regular grid, i.e., a
balanced octree with all octants at $D_{\max}$. Note that any octree
that can be generated at the same value of $D_{\max}$ can be obtained
by applying a series of local coarsening operations to the regular
grid. Since we only store the offsets it is sufficient to analyze the
distribution of the offset values for one given direction, say for a
neighbor along the $x$-axis. The expression for all other directions
are similar.

%Let us start with the root octant, i.e.,
%$d_{\max} = 0$, This is shown in the following figure,

%\tmfloat{h}{big}{figure}{\epsfig{file=}  \epsfig{file=}}{The effect of splits
%on the offsets that can be created for a regular grid.}

For $D_{\max} = 0$, there is only one octant and correspondingly the offset
is $1$. If we introduce a split in the root octant, $D_{\max}$ becomes 1, the 
offset increases by 2 for one octant. 
%This corresponds to the boundary introduced by the previous split.
On introducing further splits, the offset is going to increase for
those octants that lie on the boundaries of the original splits, and
the general expression for the maximum offset can be written as $
\tmop{offset} = 1 + \sum_{i = 1}^{D_{\max}} 2^{ d\cdot i - 1}$, for a $d$-tree. 
In addition, a number of other large offsets are produced for
intermediate split boundaries. Specifically for a regular grid at
maximum depth $D_{\max}$, we shall have $2^{d\cdot( D_{\max} - x
)}$ octants with an offset of $1 + \sum^x_{i = 1} 2^{d\cdot i - 1}$
.  As can be clearly seen from the expression, the distribution of the
offsets is geometric. With the largest number of octants having small
offsets.

For the case of general balanced octrees, we observe that
any of these can be obtained from a regular grid by a number of
coarsening operations. The only concern is whether the coarsening can
increase the offset for a given octant. The coarsening does not affect
octants that are greater than the current octant (in the Morton
order). For those which are smaller, the effect is minimal since
every coarsening operation reduces the offsets that need to be stored.

Golomb-Rice coding \cite{golomb66, rice79} is a form of entropy encoding that is
optimal for geometric distributions, that is, when small values are
vastly more common than large values. Since, the distribution of the
offsets is geometric, we expect a lot of offsets with small
values and fewer occurrences of large offsets. The Golomb coding uses 
a tunable parameter $M$ to divide an input value into two parts: 
$q$, the result of a division by $M$, and $r$, the remainder. 
In our implementation, the remainder is stored as a
\texttt{byte}, and the quotient as a \texttt{short}. On an average, we observe 
one large jump in the node indices, and therefore the amortized
cost of storing the compressed lookup table, is 8 \texttt{bytes} for
storing the remainders, 2 \texttt{bytes} for the quotient, one \texttt{byte}
for storing a flag to determine which of the 8 nodes need to use a
quotient, and one additional \texttt{byte} for storing additional element
specific flags. Storing the lookup explicitly would require 8 \texttt{ints}, 
and therefore we obtain a 3{\texttt x} compression in storing the
lookup table.

%}}}

%{{{ fem
\section{Finite Element Computation on Octrees}
\label{sec:femMatVec} 
% \label{sec:shapeFunctions}

In this section, we describe the evaluation of a \texttt{MatVec} with the global finite element `stiffness' matrix. The \texttt{MatVec} refers to
 a function that takes a vector and returns another vector, the result of applying the
  discretized PDE operator to the the input vector. The key features of our framework for performing finite element computations are listed below.
\begin{itemize}
\item We do not store hanging nodes explicitly. They do not represent independent
  degrees of freedom in a FEM solution. A method to eliminate hanging nodes in
  locally refined quadrilateral meshes and yet ensure inter-element continuity by the
use of special bilinear quadrilateral elements was presented in \cite{wwang01}. We
 extended that approach to three dimensions. If the i-th vertex of an element/octant is hanging,
 then the index corresponding to this node will point to the i-th vertex of the parent\footnote{The 2:1
 balance constraint ensures that the nodes of the parent can never be hanging.} of this element instead (Figure \ref{fig:nodes}). Thus,
 if a hanging node is shared between 2 or more elements, then in each element it might point to a different index.
\item Since we eliminate hanging nodes in the meshing stage itself, we don't need to use 
projection schemes like those used in \cite{gordon03,jones05,kim02,tu-etal06b} to enforce conformity. Hence, we don't need multiple 
tree traversals for performing each MatVec; instead, we perform a single traversal by mapping each
octant/element to one of the pre-computed element types, depending on
the configuration of hanging nodes for that element. 
\item To reduce the memory overhead, the linear octree is stored in
a compressed form that requires only one byte per octant (the level of
the octant). Even the element-to-node mappings can be compressed at a 
modest expense of uncompressing this on the fly while looping over the 
elements to perform the finite element MatVecs.
\end{itemize}

Below, we list some of the properties of the shape functions defined on octree meshes. 
\begin{itemize}
\item The shape functions are not rooted at the hanging nodes. 
\item The shape functions are tri-linear.
\item The shape functions assume a value of 1 at the vertex at which they are rooted and a value of 0 at all other non-hanging vertices in the octree.
\item The support of a shape function can spread over more than 8 elements.
\item If a vertex of an element is hanging, then the shape functions rooted at the other non-hanging vertices in that element do not vanish on this hanging node. Instead, they will vanish at the non-hanging vertex that this hanging node is mapped to. For example, in Figure \ref{fig:nodes} the shape function rooted at node $a_0$ will not vanish at nodes $a_1$, $a_2$, $a_3$, $a_4$, $a_5$ or $a_6$. It will vanish at nodes $p_1$, $p_2$, $p_3$, $p_4$, $p_5$, $p_6$ and $a_7$. It will assume a value equal to 1 at node $a_0$.
\item A shape function assumes non-zero values within an octant if and only if it is rooted at some non-hanging vertex of this octant or if some vertex of the  octant under consideration is hanging, say the $i$-th vertex, and the shape function in question is rooted at the $i$-th non-hanging vertex of the parent of this octant. Hence, for any octant there are exactly eight shape functions that do not vanish within it and their indices will be stored in the vertices of this octant.
\item The finite element matrices constructed using these shape functions are mathematically equivalent to those obtained using Projection schemes such as in \cite{kim02, tu05, tu-etal06b}.
\end{itemize}

To implement finite element MatVecs using these shape functions, we need to enumerate all the permissible hanging configurations for an octant. The following properties of 2:1 balanced linear octrees helps reduce the total number of permissible hanging configurations. Figure \ref{fig:nodes} illustrates these properties.
\begin{itemize}
\item Every octant has at least 2 non-hanging nodes and they are:
\begin{itemize}
\item The node that is common to both this octant and its parent.
\item The node that is common to this octant and all its siblings.
\end{itemize}
\item An octant can have a face hanging node only if the remaining nodes on that face are one of the following:
\begin{itemize}
\item Edge hanging nodes.
\item The node that is common to both this octant and its parent.
\end{itemize}
\end{itemize}

% \begin{figure}
%   \begin{center} 
%   \input{TIKZ_Figures/hanging} 
%   \caption{\label{fig:hanging} Illustration of nodal-connectivities required to perform conforming
%   FEM calculations using a single tree traversal. Every octant has at least 2 non-hanging nodes, one of which is
%   shared with the parent and the other is shared amongst all the
%   siblings. The octant shown in blue $(a)$ is a child 0, since it shares its zero node
%   $(a_0)$ with its parent $(p)$. It shares node $a_7$ with its siblings. All other nodes, if hanging, 
%   point to the corresponding node of the parent octant instead. Nodes, $a_3, a_5, a_6$ are face
%   hanging nodes and point to $p_3, p_5, p_6$, respectively. Similarly
%   $a_1, a_2, a_4$ are edge hanging nodes and point to $p_1, p_2,
%   p_4$. All the nodes in this illustration are labelled in the Morton ordering. } 
%   \end {center}
% \end{figure}

The 'child number' of an octant is another important property that is used to handle hanging nodes in this framework. This is used both while implementing the MatVec for the finite element matrices and for implementing the inter-grid transfer operations. An octant's position with respect to its parent is identified by specifying the node that it shares with its parent. If an octant shares its k-th node with its parent, then it is said to have a 'child number' equal to k. For convenience, we use the Morton ordering to number the nodes of an octant. Thus, sorting the children of an octant in the Morton order is equivalent to sorting the children according to their child numbers. The child number of an octant is a function of the coordinates of its anchor and its level in the tree. Algorithm \ref{alg:childNumber} is used to compute the child number of an octant.

\begin{table*}
\centering
\rule{\textwidth}{0.01mm}
\begin{algorithm}
{\textsc{Finding the child number of an octant}} \\
\rule{\textwidth}{0.01mm}
\flushleft
\tt{\bf{Input:~} The anchor (x,y,z) and level (d) of the octant and the maximum permissible depth of the tree ($D_{max}$).} \\
\tt{\bf{Output:} $c,$ the child number of the octant.}\\
~\\
1. $l \leftarrow 2^{\left( D_{max} - d \right)}$\\
2. $l_p \leftarrow 2^{\left( D_{max} - d + 1 \right)}$ \\
3. $(i,j,k) \leftarrow (x,y,z) \bmod l_p$ \\
4. $(i,j,k) \leftarrow (i,j,k) / l$ \\
5. $c \leftarrow \left(4k + 2j + i\right)$\\
\label{alg:childNumber}
\end{algorithm}
\rule{\textwidth}{0.01mm}
\end{table*}

In each of the 8 child number based configurations (Figures \ref{fig:child0} - \ref{fig:child7}), $v_0$ is the node that this element shares with its parent and $v_7$ is the node that this element shares with all its siblings. For an element with child number $k$, $v_0$ will be the k-th node and $v_7$ will be the $(7-k)$-th node. $v_0$ and $v_7$ can never be hanging. If $v_3$, $v_5$ or $v_6$ are hanging, they will be face-hanging and not edge-hanging. If $v_3$ is hanging, then $v_1$ and $v_2$ must be edge-hanging. If $v_5$ is hanging, then $v_1$ and $v_4$ must be edge-hanging. If $v_6$ is hanging, then $v_2$ and $v_4$ must be edge-hanging. After factoring in these constraints, there are only 18 potential hanging node configurations for each of the 8 'child number' configurations (Table \ref{tab:hnTypes}).

\begin{figure}
  \begin{center} 
  \subfloat[Child 0] {  
  \includegraphics{TIKZ_Figures/child0} 
  } 
  \subfloat[Child 1] {
   \label{fig:child1}
  \includegraphics{TIKZ_Figures/child1} 
  }  
  \subfloat[Child 2] {
   \label{fig:child2}
  \includegraphics{TIKZ_Figures/child2} 
  }  \\
  \subfloat[Child 3 ] {
   \label{fig:child3}
  \includegraphics{TIKZ_Figures/child3} 
  }      
  \subfloat[Child 4] {
   \label{fig:child4}
  \includegraphics{TIKZ_Figures/child4} 
  } 
  \subfloat[Child 5] {
   \label{fig:child5}
  \includegraphics{TIKZ_Figures/child5} 
  } \\
  \subfloat[Child 6]{
   \label{fig:child6}
  \includegraphics{TIKZ_Figures/child6} 
  }
  \subfloat[Child 7]{
   \label{fig:child7}
  \includegraphics{TIKZ_Figures/child7} 
  }
  \caption{Illustration of the 8 child number configurations. In each of these configurations, $v_0$ is the node that this element shares with its parent and $v_7$ is the node that this element shares with all its siblings. For an element with child number $k$, $v_0$ will be the k-th node and $v_7$ will be the $(7-k)$-th node. $v_0$ and $v_7$ can never be hanging. If $v_3$, $v_5$ or $v_6$ are hanging, they will be face-hanging and not edge-hanging. If $v_3$ is hanging, then $v_1$ and $v_2$ must be edge-hanging. If $v_5$ is hanging, then $v_1$ and $v_4$ must be edge-hanging. If $v_6$ is hanging, then $v_2$ and $v_4$ must be edge-hanging. All the octants shown are oriented in the global coordinate system.} 
  \end {center}
\end{figure}
    
\begin{table}
	\centering
		\begin{tabular}{|c|l|}\hline
		\multicolumn{1}{|c}{Element Type}  & \multicolumn{1}{|c|}{Corresponding Hanging nodes}\\\hline
		    & \\
		 1     & None\\
		 &  \\
     2    & $v_2$\\
     &  \\
     3   & $v_1$\\
     &  \\
     4    & $v_1$ and $v_2$\\
     &  \\
     5    & $v_4$\\
     &  \\
     6    & $v_4$ and $v_2$\\
     &  \\
     7    & $v_4$ and $v_1$\\
     &  \\
     8    & $v_4$, $v_1$ and $v_2$\\
     &  \\
     9    & $v_3$, $v_1$ and $v_2$\\
     &  \\
     10   & $v_3$, $v_4$, $v_1$ and $v_2$\\
     &  \\
     11    & $v_6$, $v_4$ and $v_2$\\
     &  \\
     12    & $v_6$, $v_4$, $v_1$ and $v_2$\\
     &  \\
     13    & $v_6$, $v_3$, $v_4$, $v_1$ and $v_2$ \\
     &  \\
     14     & $v_5$, $v_4$ and $v_1$\\
     &  \\
     15      & $v_5$, $v_4$, $v_1$ and $v_2$\\
     &  \\
     16  & $v_5$, $v_3$, $v_4$, $v_1$ and $v_2$\\
     &  \\
     17    & $v_5$, $v_6$, $v_4$, $v_1$ and $v_2$ \\
     &  \\
     18   & $v_5$, $v_6$, $v_3$, $v_4$, $v_1$ and $v_2$ \\
     &  \\\hline     
		\end{tabular}
\caption{The list of permissible hanging node configurations for any octant.}
\label{tab:hnTypes}
\end{table}

\subsection{Overlapping communication with computation}
\label{sec:overlapCommCompFE}
Every octant is owned by a single processor. However, the values of unknowns associated with octants on inter-processor boundaries need to be shared among several processors. We keep multiple copies of the information related to these octants and we term them `ghost' octants. In our implementation of the finite element MatVec, each processor iterates over all the octants it owns and also loops over a layer of ghost octants that contribute to the nodes it owns. Within the loop, each octant is mapped to one of the above described hanging configurations. This is used to select the appropriate element stencil from a list of pre-computed stencils. We then use the selected stencil in a standard element based assembly technique. Although the processors need to read ghost values from other processors they only need to write data back to the nodes they own and do not need to write to ghost nodes. Thus, there is only one communication step within each MatVec. We even overlap this communication with useful computation. We first loop over the elements in the interior of the processor domain since these elements do not share any nodes with the neighbouring processors. While we perform this computation, we communicate the ghost values in the background. At the end of the first loop, we use the ghost values that were recieved from the other processors and loop over the remaining elements.



%These eight types are identified by the number of face hanging and edge hanging nodes in the element under consideration: (0,0),(0,1),(0,2),(0,3),(1,2),(1,3),(2,3), and (3,3) where the first number is the number of face hanging nodes and the second number is the the number of edge hanging nodes.

The \texttt{MatVec} can be evaluated either by looping through the elements or by looping through the
 nodes. Due to the simplicity of the elemental loop, we have
chosen to support elemental loops over nodal loops. We first loop over the 
elements in the interior of the processor domain since these elements do not share 
any nodes with the neighboring processors. While we perform this computation we communicate 
the values of the ghost nodes in the background. At the end of the first loop, we use the ghost 
values from other processors and loop over the remaining elements.

%In the case of an elemental loop, we need to store the information pertaining
%to the eight nodes defining this hexahedral element. In the case of a
%nodal loop, we need to store the information about 26
%neighboring nodes.
%}}}

\section{Summary of Octree Algorithms}
\label{sec:octsummary}
The sequence of steps involved in solving a PDE on octrees is summarized below:

\begin{enumerate}
\item A `sufficiently' fine\footnote{Here the term sufficiently is used to mean that the discretization error introduced is acceptable.} 2:1 balanced complete linear octree is constructed using the algorithms described in Section \ref{sec:methods}
\item The octree at is meshed using the algorithm described in Section \ref{sec:oda}. 
\item The discrete system of equations is then solved using the conjugate gradient algorithm.
\end{enumerate}

Table \ref{tab:complexity} gives the parallel time complexity of the various algorithmic components as a function of the problem size, $N$, and the number of processors, $n_p$.

\begin{table}   
  \begin{center}
  \begin{tabular}{|l|l|} \hline 
  {Algorithmic} &  {Parallel}\\
  {Component} & {Time Complexity} \\ \hline
  {Octree Construction}  & $\mathcal{O}(\frac{N}{n_p}\log(\frac{N}{n_p}) + n_p \log n_p)$\\
  &\\
  {2:1 Balancing}  & $\mathcal{O}(\frac{N}{n_p}\log(\frac{N}{n_p}) + n_p \log n_p)$\\
  &\\
  {Block Partition} & $\mathcal{O}(\frac{N}{n_p} + n_p)$\\
  &\\
  {Meshing}  & $\mathcal{O}(\frac{N}{n_p}\log(\frac{N}{n_p}) + n_p \log n_p)$\\
  &\\
  {Matvecs}  & $\mathcal{O}(\frac{N}{n_p})$\\ \hline
  \end{tabular}
\end{center}  
\caption{Complexity estimates for the various algorithmic components assuming a Hypercube network topology with $\theta(n_p)$ bandwidth. $N$ is the size of the linear octree and $n_p$ is the number of processors.}
\label{tab:complexity}
 \end{table}


